Sex & statistical power - simulations

Code

library(ggplot2)
library(pander)
library(emmeans)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ lubridate 1.9.2     ✔ tibble    3.1.8
✔ purrr     1.0.1     ✔ tidyr     1.3.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

library(here)

here() starts at /Users/szymek/Library/CloudStorage/Dropbox/#SCIENCE/R/230713_power_sex_sim

Code

library(ggpubr)

Analysis preps (colours, options, auxiliary functions)

Code

cbb_palette <- c(
  "#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
  "#0072B2", "#D55E00", "#CC79A7"
)

## Avoid table wrapping
options(width = 800)
set.seed(999)

source(here("PhillipsKarpEtal_data/funs.r"))

Below results repeat simulations from the paper of Phillips et al. (2018) PLoS Biol 21(6): e3002129. Original results simplified by removing the post-hoc and the (pooled) t-test methods of analysing data. Generation of simulated data was modified to use a formal linear model defined as:

\[\mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e}\]

with \(\mathbf{b}=(\mu,\beta_{Treat},\beta_{Sex},\beta_{Treat:Sex})\) and \(\mathbf{e} \sim \mathcal{N}(\mathbf{0}, \mathbf{s})\), with \(\mathbf{s}=\sigma\sqrt{\mathbf{1}+\mathbf{X}_{Sex}\cdot\alpha}\) (\(\alpha\) quantifies the amount of heteroscedasticity in sex variances, with SD \(\sigma_{m}=\sigma_{f}\sqrt{1+\alpha}\), or male variance being \(\alpha\cdot100\%\) greater than female). Regression coefficients (\(\mathbf{b}\)) are compatible with contrasts (sex/treatment effects) used in Phillips et al. (2023).

1 Scenarios from Fig. 3

1.1 From original paper

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none <- all_effs_sex_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

calc_power_sex_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1

all_effs_sex_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small <- all_effs_sex_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none"))

calc_power_sex_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1

all_effs_sex_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large <- all_effs_sex_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none"))


calc_power_sex_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_none, calc_power_sex_small,
                               calc_power_sex_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
             group = interaction(sex_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(sex_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sex_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sex effect size") + ylim(0, 1)

1.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none_lowh <- all_effs_sex_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small_lowh <- all_effs_sex_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low"))

calc_power_sex_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_none_largeh <- all_effs_sex_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_small_largeh <- all_effs_sex_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

It should be of note that evolutionary biology has seen heteroscedasticities exceeding 100% in variance ration between sexes. Combined results below show the impact of heteroscedasticity.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_small,
                               calc_power_sex_large,
                               calc_power_sex_small_lowh,
                               calc_power_sex_large_lowh,
                               calc_power_sex_small_largeh,
                               calc_power_sex_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(sex_eff_size = factor(sex_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
             group = interaction(sex_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

plotA <- all_sex_effs_sex_only %>%
  filter(sex_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

plotA

2 Scenarios from Fig. 4 - interaction of magnitude

2.1 From original paper

Below simulations were only slightly modified. We assumed a moderate (+0.5) overall effect of sex and a varying interaction increasing the magnitude of sex difference in the treated group (+0, +0.5, +1).

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_ix_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1

all_effs_sex_ix_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1

all_effs_sex_ix_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("large"))


calc_power_sex_ix_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
                               calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(ix_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = ix_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)

2.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_ix_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_ix_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Power consequences of heteroscedasticities with varying (magnitude-only) interaction.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
                               calc_power_sex_ix_large,
                               calc_power_sex_ix_small_lowh,
                               calc_power_sex_ix_large_lowh,
                               calc_power_sex_ix_small_largeh,
                               calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

plotB <- all_sex_effs_sex_only %>%
  filter(ix_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

plotB

3 Scenarios from Fig. 5 - interaction of direction

3.1 From original paper

Below simulations were only slightly modified. We assumed a strong (-1) overall effect of sex and a varying interaction resulting of crossing of sex effects on the treatment effect (0, -1, -2).

Assuming homoscedasticity of sex variances

No sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0

all_effs_sex_ix_none <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1

all_effs_sex_ix_small <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1

all_effs_sex_ix_large <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("none")) %>%
    mutate(ix_eff_size = rep("large"))


calc_power_sex_ix_large %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Plots combining all scenarios from Fig. 1

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
                               calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, method))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = method)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Stat method")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  # filter(!Effect %in% c("T-Test Treatment Eff")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
             group = interaction(ix_eff_size, Effect))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = ix_eff_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)

3.2 Heteroscedastic sex variances

We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.

Low heteroscedasticity

Code

alpha <- 1

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_lowh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha

all_effs_sex_ix_small_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
  mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("low")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_lowh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity

Code

alpha <- 2

Small sex eff, no ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_ix_none_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
    filter(!Effect == "Residuals    ") %>%
    group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
    summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
    mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("none"))

calc_power_sex_ix_none_largeh %>%
    ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
    geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
    geom_line(aes(group = Effect)) +
    theme_classic() +
    ylab("Power") + ylim(0, 1)

Small sex eff, small ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha

all_effs_sex_ix_small_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("small")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("small"))

calc_power_sex_ix_small_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Small sex eff, large ix

Code

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(heterosc = rep("large")) %>%
    mutate(ix_eff_size = rep("large"))

calc_power_sex_ix_large_largeh %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Power consequences of heteroscedasticities with varying (direction-only) interaction. Power anomalies indicate pretty low inferential abilities of such models at low sample sizes.

Code

all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
                               calc_power_sex_ix_large,
                               calc_power_sex_ix_small_lowh,
                               calc_power_sex_ix_large_lowh,
                               calc_power_sex_ix_small_largeh,
                               calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

all_sex_effs_sex_only %>%
  filter(Effect %in% c("Treatment")) %>%
  mutate(ix_eff_size = factor(ix_eff_size,
                               levels = c("small", "large"))) %>%
  mutate(method = recode(method,
                         "ANOVA" = "Factorial")) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
             group = interaction(ix_eff_size, heterosc))) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity")  + ylim(0, 1)

Code

all_sex_effs_sex_only %>%
  filter(ix_eff_size %in% c("large")) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  mutate(heterosc = factor(heterosc,
                               levels = c("none", "low", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = heterosc)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)

4 Impact of sample size on power

4.1 Large sex eff, no ix

Below simulation code generates power estimates for the following sample size scenarios: 5M+5F in each treatment group (original scenario, totql sample 20 individuals), 10M+10F (N = 40), 50M+50F (N = 200).

First - low heteroscedasticity:

5M+5F

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_5 <- all_effs_sex_large_lowh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
  mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_large_lowh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_10 <- all_effs_sex_large_lowh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_large_lowh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_lowh_50 <- all_effs_sex_large_lowh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_large_lowh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity:

5M+5F

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_5 <- all_effs_sex_large_largeh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_large_largeh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_10 <- all_effs_sex_large_largeh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_large_largeh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_large_largeh_50 <- all_effs_sex_large_largeh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("none")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_large_largeh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

4.2 Large sex effect, crossed reaction norms (interaction)

Identical sample size scenarios.

First - low heteroscedasticity:

5M+5F

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_lowh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_5 <- all_effs_sex_ix_large_lowh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_ix_large_lowh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_10 <- all_effs_sex_ix_large_lowh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_ix_large_lowh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_lowh_50 <- all_effs_sex_ix_large_lowh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("low"))

calc_power_sex_ix_large_lowh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

Large heteroscedasticity:

5M+5F

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha

all_effs_sex_ix_large_largeh_5 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_5 <- all_effs_sex_ix_large_largeh_5 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("5M+5F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_5 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

10M+10F

Code

no_per_gp <- 10

all_effs_sex_ix_large_largeh_10 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_10 <- all_effs_sex_ix_large_largeh_10 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("10M+10F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_10 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

50M+50F

Code

no_per_gp <- 50

all_effs_sex_ix_large_largeh_50 <- rbind(
  p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
    )
)

calc_power_sex_ix_large_largeh_50 <- all_effs_sex_ix_large_largeh_50 %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
    mutate(ix_eff_size = rep("large")) %>%
    mutate(sample_size = rep("50M+50F")) %>%
    mutate(heterosc = rep("large"))

calc_power_sex_ix_large_largeh_50 %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(group = Effect)) +
  theme_classic() +
  ylab("Power") + ylim(0, 1)

4.3 Summaries for sample size scenarios

No interaction, low (left) and large (right) heteroscedasticity:

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_lowh_5,
    calc_power_sex_large_lowh_10,
    calc_power_sex_large_lowh_50,
    calc_power_sex_large_largeh_5,
    calc_power_sex_large_largeh_10,
    calc_power_sex_large_largeh_50
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("low")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)


plot2 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("large")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)

ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

Interaction of effect direction present, low (left) and high (right) heteroscedasticity:

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_ix_large_lowh_5,
    calc_power_sex_ix_large_lowh_10,
    calc_power_sex_ix_large_lowh_50,
    calc_power_sex_ix_large_largeh_5,
    calc_power_sex_ix_large_largeh_10,
    calc_power_sex_ix_large_largeh_50
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("low")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)


plot2 <- all_sex_effs_sex_only %>%
  filter(heterosc %in% c("large")) %>%
  mutate(sample_size = factor(sample_size,
                               levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
  geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line(aes(linetype = sample_size)) +
  theme_classic() +
  ylab("Statistical power") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) +
  scale_linetype_discrete(name = "Sample size") + ylim(0, 1)

ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

5 Simulating sample size

5.1 Simulations assuming target power = 80%

Scenario 1: strong sex effect, none to large heteroscedasticity. Assumed power 80% for treatment effect.

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot1

Code

plot1.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 70)

plot1.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot2 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot2

Code

plot2.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 25)

plot2.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

Code

alpha <- 0

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_noh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("none"))

Code

alpha <- 1

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_lowh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("low"))

Code

alpha <- 2

no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha

all_effs_sex_large_largeh <- rbind(
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    ),
    p_value_interaction_crossed(
        n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
        variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
        male_es = male_es, female_es = female_es,
        treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
        fix_power = T
    )
)

calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
  filter(!Effect == "Residuals    ") %>%
  group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
  summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
  mutate(sex_eff_size = rep("large")) %>%
  mutate(ix_eff_size = rep("none")) %>%
    mutate(heterosc = rep("large"))

Code

all_sex_effs_sex_only <- rbind(
    calc_power_sex_large_noh,
    calc_power_sex_large_lowh,
    calc_power_sex_large_largeh
)

all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))

# plot1 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)


# plot2 <- all_sex_effs_sex_only %>%
#   filter(Effect %in% c("Treatment")) %>%
#   filter(heterosc %in% c("large")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(method = recode(method,
#                          "ANOVA" = "Factorial")) %>%
#   ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
#   geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
#   geom_line() +
#   theme_classic() +
#   ylab("Statistical power") +
#   xlab("Main treatment effect") +
#   scale_colour_manual(name = "Sample size", values = cbb_palette) +
#   # scale_linetype_discrete(name = "Heteroscedasticity")  +
#   ylim(0, 1)

# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")

plot3 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
  theme(text = element_text(size=16))

plot3

Code

plot3.1 <- all_sex_effs_sex_only %>%
#   filter(heterosc %in% c("low")) %>%
#   mutate(sample_size = factor(sample_size,
#                                levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
#   mutate(sex_eff_size = factor(sex_eff_size,
#                                levels = c("none", "small", "large"))) %>%
  ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
  # geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
  geom_line() +
  theme_bw() +
  ylab("Sample size") +
  xlab("Main treatment effect") +
  scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
  xlim(0.3, 0.7) + ylim(5, 2200)

plot3.1

Warning: Removed 36 rows containing missing values (`geom_line()`).

5.2 Combined plots

(In all below examples the y axis (sample size) indicates the number of individuals per group and sex - hence the total sample size is 4 times the number indicated on the y axis.)

Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for full range of effect sizes and on log-scale for sample size:

Code

comb1 <- ggarrange(plot1, plot2, plot3, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))

comb1

Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for moderate-sized effect sizes using identity scale for sample size:

Code

comb2 <- ggarrange(plot1.1, plot2.1, plot3.1, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))

Warning: Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).
Removed 36 rows containing missing values (`geom_line()`).

Code

comb2

Combined plots demonstrating the impact of heteroscedasticity on power for a no-interaction case (A) and interaction-of-magnitude case (B; stronger effect in the more variable sex):

Code

comb3 <- ggarrange(plotA, plotB, ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom", labels = c("A", "B"))

comb3